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We discuss the initial-boundary value problem for the Baumgarte-Shapiro-Shibata-Nakamura evo- 
lution system of Einstein's field equations which has been used extensively in numerical simulations 
of binary black holes and neutron stars. We specify nine boundary conditions for this system with 
the following properties: (i) they impose the momentum constraint at the boundary, which is shown 
to preserve all the constraints throughout evolution, (ii) they approximately control the incoming 
gravitational degrees of freedom by specifying the Weyl scalar at the boundary, (iii) they control 
the gauge freedom by requiring a Neumann boundary condition for the lapse, by setting the normal 
CNj ■ component of the shift to zero, and by imposing a Sommerfeld-like condition on the tangential com- 

ponents of the shift, (iv) they are shown to yield a well-posed problem in the limit of weak gravity. 
' ' Possible numerical applications of our results are also discussed briefly. 

PACS numbers: 04.20.Ex, 04.25. D-, 04.20.-q 



I. INTRODUCTION 
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Many of the boundaries introduced in the modeling of physical problems are not part of the system. In the 
numerical Cauchy evolution of a gravitational system, for instance, one is faced with the finiteness of the numerical 
| grid which usually leads to an evolution scheme on a spacetime domain with artificial initial and outer boundary 
surfaces. Initial conditions are chosen such that they represent as accurately as possible the physical system at one 
' moment of time. The outer boundary surface in turn needs to be transparent to the system in the sense that the 
\Q , conditions imposed on this surface mimic as well as possible the evolution on an infinite, asymptotically flat domain. 
This leads to the construction of absorbing outer boundary conditions which can be defined by the requirement 
that they yield a well-posed initial-boundary value problem (IBVP) and minimize spurious reflections from the outer 
f^i • boundary. Furthermore, the boundary conditions should be chosen such that they preserve the constraints throughout 
t— I | evolution. 

There has already been a large amount of work on outer boundary conditions in General Relativity, see [l| for a 
recent review. Absorbing boundary conditions which preserve the constraints and lead to a well-posed IBVP have first 
, been constructed for a tetrad formulation Q and more recently also for the harmonic formulation of Einstein's 
vacuum field equations. In particular, error estimates for spurious reflections and hig her-order absorbing boundary 
rS | conditions have been given in 0, H| and implemented and tested numerically in 0, Il0j |. 

However, these results are not yet applicable to the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation 
[ill [l2j of Einstein's field equations on which many numerical codes for simulating binary black holes or neutron 
stars are based. So far, such codes use ad-hoc "radiative boundary condition" (see, for example, Ref. [13), which 
are imposed on all the evolution variables. Although these conditions seem very easy to implement numerically, they 
are unlikely to yield a well-posed Cauchy evolution or to preserve the constraints. This results in a lack of either 
accuracy or efficiency since the outer boundary has to be pushed very far away in the wave zone in such a way that 
the boundary surface is causally disconnected from the region where physics is extracted. 

In this article we construct absorbing boundary conditions for the BSSN system which preserve the constraints and 
discuss the well-posedness of the resulting initial-boundary value problem. Previous work on specifying boundary 
conditions for the BSSN system has been carried out in [14], where the BSSN system for a rather general class 
of gauge conditions on the lapse and a fixed shift tangent to the boundary is shown to be reducible to a first- 
order symmetric hyperbolic system (FOSH). Then, boundary conditions are given to the six incoming characteristic 
fields, resulting in a well-posed IBVP. However, the conditions presented in [l4| do not preserve the constraints. 
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As a consequence, constraint-violating modes may be generated at the boundary surface. The work in [l6j on the 
other hand, formulates constraint-preserving boundary conditions for the BSSN system but the well-poscdncss of the 
resulting IBVP in the absorbing case is not established. Furthermore, numerical simulations of binary black hole 
mergers require more general gauge conditions on the shift than the ones considered for the IBVPs formulated in 

Mm 

Here, we present a new set of boundary conditions for the BSSN system with the "hyperbolic i^-driver" condition for 
the lapse and the "hyperbolic Gamma-driver" condition for the shift vector [ItJ which is currently used in numerical 
simulations of binary black holes. The main properties of our boundary conditions are: (i) they preserve the constraints 
throughout evolution, (ii) they control the Weyl scalar ^ at the boundary, a condition that yields small spurious 
reflections of gravitational radiation 0, |1[ , (iii) they provide conditions on lapse and shift in the form of a Neumann 
condition for the lapse, a Dirichlet condition for the normal component of the shift and a Sommerfeld-like boundary 
condition on the tangential components of the shift, (iv) they are shown to yield a well-posed IBVP in the case of 
weak gravitational fields, where the equations may be linearized about Minkowski spacetime. If the outer boundary 
is placed in the weak field zone, this should be a reasonable assumption, and in view of the results in Ref. [3l about 
the well-posedness of the Cauchy problem without boundaries one might except that our boundary conditions yield 
a well-posed IBVP in the nonlinear case as well. 

Our boundary conditions are summarized in Sec. |H] along with the BSSN evolution equations and constraints, and 
in the rest of the article the derivation of these boundary conditions is presented. The derivation is based on the 
analysis of small amplitude, high-frequency perturbations of the system which gives rise to a linear system on the 
half space x > with planar boundary at x = 0. Boundary conditions are first constructed for the linearized system 
and then extrapolated to the nonlinear case. The linear system is studied in Sec. IIIII using the methods described 
in [18, fl9j where boundary conditions are constructed in three steps. In a first step, the propagation of the gauge- 
invariant guantities, consisting of the constraint variables and the linearized Weyl curvature is analyzed. It is shown 
that this system yields a FOSH evolution system. We specify boundary conditions for this system which yield a 
well-posed IBVP, guarantee the propagation of the constraint fields and also control the Weyl scalar ^o- The second 
step consists in controlling the gauge degrees of freedom. The gauge functions, the lapse a, and the shift vector /3 l are 
free to be chosen as better fits the problem. That freedom, however, is diminished by several practical requirements. 
The hyperbolic K-drrver and Gamma-driver mentioned above are used in most BSSN formulations. We show that 
this system can be cast into FOSH form, and specify well-posed boundary conditions. In particular, the boundary 
conditions require that the normal component of the shift is zero, a condition that is important in order to ensure that 
the number of ingoing characteristic fields is constant throughout evolution, as explained in the next section. The final 
step consists in reconstructing the BSSN variables and to show that the initial data and our boundary conditions give 
rise to unique solutions of the IBVP. The residual gauge freedom in the linearized IBVP is also analyzed in Sec. IIIII 
and it is shown to be parametrized by two functions and three vector fields on the initial slice and one transverse 
vector field on the boundary. 

Next, in Sec. IIVI we show that our boundary conditions do, in fact, preserve the constraints in the nonlinear 
case. This means that solutions of the nonlinear BSSN system satisfying our boundary conditions and the initial 
constraints satisfy the constraints everywhere on the spacetime domain. We present a brief discussion of a possible 
numerical implementation in Sec. |V] and draw conclusions in Sec. I VII Since our methods are based on the theory of 
FOSH systems with maximal dissipative boundary conditions, we summarize the relevant results of this theory in an 
appendix. 

We finish our introduction by mentioning that there are alternative methods to cope with the problem of representing 
an infinite domain on a finite numerical grid. Examples are the Cauchy-characteristic matching method [20] and the 
hyperboloidal initial value problem where a compactification along hyperboloidal surfaces in a scri-fixing gauge allows 
to describe the gravitational waveform at null infinity in an unambiguous way, and no outer boundary conditions are 
needed. We refer the reader to Refs. [2l| - [26j for recent developments and the original references on the hyperboloidal 
initial value problem. While still under development, it constitutes a very promising approach for the computation 
of gravitational waves emitted by isolated systems. 

II. BOUNDARY CONDITIONS FOR THE BSSN SYSTEM 

The evolution equations we are referring to in our work are the BSSN evolution equations with the hyperbolic 
if-driver and Gamma-driver gauge conditions as described in 17] with the exception of the advection terms. Using 
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the notation in [14J this system reads 
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where we have decomposed the three metric and the extrinsic curvature according to 
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7y having unit determinant and Aij being trace- free, and where we have introduced the operator 8q — dt — P°dj. 
Here, all quantities with a tilde refer to the conformal three metric 7y, and the latter is used in order to raise and lower 
their indices. In particular, Di and T k ^ refer to the covariant derivative and the Christoffel symbols, respectively, with 
respect to 7^ . The expression [...] denotes the trace-less part (with respect to the metric 7y) of the expression 
inside the parentheses, and 
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R% = -2D i D j( f>-2^ ij D k Dk4> + 4:D i (l>D j 4>-Aj ij D k ( l)D k (f>. 

The parameter m, which was introduced in .27], controls how the momentum constraint is added to the evolution 
equations for the variable T l . Also, Ko is a smooth function of its argument, /, G and H are strictly positive and 
smooth functions of their arguments, and rf is a vector-valued function depending smoothly on its arguments. The 
choice 



m=l, f(a,<P,x fl ) = ~ , tf (O = 0, 

a 



r) i (B j ,a,x> J ')=r)B i , 



with rj a positive constant corresponds to the evolution system used in recent black hole simulations based on 1 + log 
slicing and the moving puncture technique (see, for instance, Ref. (28|) or the turducken approach [29|, |3(|. The 
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TABLE I: The nine boundary conditions for the BSSN variables. Here, n is the unit outward normal to the boundary which 
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tangential vectors and trace-less tangential tensors, respectively. We are also defining the quantities Sij = Rij + R^ij + 
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which determine the electric and magnetic parts 



of the Weyl tensor through Eij — Sij — |7ij7 fei l £ki and Bij — Bki(i£j) ki , respectively, where Ekij denotes the volume form with 
respect to the three metric 7^. Finally, Gij is a given function on the boundary which determines the value of the Weyl scalar 
$0 at the boundary. The precise relation between Gij and is the following: if N = a~ (dt — 0'di) denotes the future-directed 
unit normal to the time slices, we may construct a Newman-Penrose null tetrad {l,k,m,fh} at the boundary by defining the 
real null vectors / := 2~ 1 ' /2 (7V + n), k := 2~ 1//2 (7V — n), and choosing m to be a unit complex null vector orthogonal to I and 
k. Then, \&o = (Eki — iBki)m k m l — Gkim k m l . For typical applications involving the modeling of isolated systems one may set 
dj to zero or freeze its value to the one computed from the initial data. 
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Boundary Conditions 



Gauge condition on the lapse 
Gauge condition on the normal component of the shift 

Gauge condition on the tangential components of the shift 

Constraint preserving condition 
^0 specifying condition 
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source terms S, Sij and S l are denned in terms of the four Ricci tensor, , and the constraint variables 1 
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The Einstein field equations are equivalent to the evolution equations 
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constant and where p, ji and sy denote, respectively, the en ergy density, the momentum density and the stress tensor 
of the matter fields as measured by Eulerian observers, see [131 ] for a definition. 

Considering those readers who already want to see the end of the story, we directly present our main result, and 
describe their derivation in detail below. Our nine boundary conditions for the BSSN system with m = 1, consistent 
with the constraint equations are given in Table [I] Let us explain first why there should be nine boundary conditions. 
In [3] the BSSN system was analyzed in terms of the eigenfields, and it was shown that under certain restrictions on 
m and the functions /, G and H the Eqs. (JTHHJ) yield a strongly hyperbolic system giving rise to a well-posed time 
evolution. For m = 1 the restrictions reduce to 



AGH f 3/. 

The characteristic speeds with respect to the time evolution vector field dt are then 
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1 Notice that these constraint variables are related to variables H, Mi and defined in [2E 
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where n refers to the unit outward normal to the boundary tangent to the time foliation. This result is based on a 
first-order pseudo-differential reduction of the equations psH |32| which does not introduce any artificial constraints. 
Therefore, these speeds are intrinsic to the BSSN formulation. If <9 t is tangent to the boundary, which is usually 
assumed in order to avoid the boundary from moving through the computational domain, the sign of these speeds 
determines the number of incoming characteristic fields and hence also the number of boundary conditions that 
must be specified. Namely, the number of boundary conditions is equal to the number of eigenfields with positive 
speed. Because of the presence of the modes with speed /3™ it is the sign of the normal component of the shift which 
determines the number of incoming fields at the boundary. The simplest way of controlling the number of boundary 
conditions throughout evolution is to demand that /3™ be zero at the boundary at all times. Otherwise, in order to 
gain control on the in- and out-going eigenfields, a careful study on the behavior of the shift vector would have to be 
implemented at each time step. Therefore, in our study, we take the simplest approach and set /3™ = as one of our 
boundary conditions. The analysis in [lij then reveals that there are precisely nine incoming eigenfields and thus, 
nine conditions have to be imposed at the boundary. If more conditions are given, the system will be overdetermined 
and there will not be solutions; if less conditions are given, the solutions will not be unique. 

As our analysis in the weak field regime shows, the nine boundary conditions are distributed as follows: there are 
four conditions that must be imposed for the gauge functions, namely the lapse and shift. One of these conditions 
sets the normal component of the shift to zero, as explained above. Geometrically, this implies that the boundary 
surface is orthogonal to the time foliation. Then, there is a Neumann boundary condition for the lapse. The other 
two gauge boundary conditions are Sommerfeld-like boundary conditions involving the tangential components of 
the shift and the tangential derivatives of the lapse. These conditions arise from the analysis of the characteristic 
structure of the gauge sector, and it would have been hard to guess them. As described in the next section, an 
alternative to these conditions is to directly specify the tangential components of the shift vector at the boundary. 
Next, there are three boundary conditions coming from the requirement, which is quite natural, of the momentum 
constraint being satisfied at the boundary. As we show in Sec. IIV1 these conditions imply constraint preservation 
without overdetermining the system. This means that the time evolution of initial data satisfying the constraints 
automatically satisfies the constraints on all time slices, and that small initial violations of the constraints which are 
usually present in numerical applications yield solutions to the evolution system where the growth of the constraint 
violations is controlled. In particular, we note that the Hamiltonian constraint must not be imposed additionally at 
the boundary, otherwise the resulting evolution system is overdetermined when small violations of the constraints 
are present. Finally, the last two boundary conditions are related to the actual gravitational degrees of freedom in 
the following sense: we are thinking about a problem which is localized and isolated and where the outer boundary 
is placed in the weak field and wave zone. In this case, the peeling behavior of the Weyl components [HI can be 
considered to be valid, and the complex Weyl scalar can be interpreted as approximately describing the incoming 
gravitational radiation. Therefore, it is natural to set V^o to zero or to its value computed from the initial data in 
order to minimize incoming gravitational radiation. While this condition only makes precise sense at null infinity, it 
has been successfully numerically implemented and tested for truncated domains with artificial boundaries, see for 
example @. Estimates on the amount of spurious reflections introduced by this condition have also been derived in 

The boundary conditions in Table U are extrapolated from the boundary conditions constructed in the next section 
for the linearized system. Some ideas about their numerical implementations are discussed in Sec. |V] 



III. ANALYSIS OF THE LINEARIZED SYSTEM 



For the following, we restrict our analysis to the case m = 1 since this is the common choice in numerical relativity 
applications. Furthermore, we consider only small amplitude, high-frequency perturbations of smooth solutions, which 
intuitively, constitute the relevant limit for analyzing the continuous dependence of the solution on the initial data [34j . 
In this limit only the principal part of the equations matters and the coefficient appearing in front of the derivative 
operators can be frozen to their value at an arbitrary point p. Therefore, the study simplifies to a linear evolution 
problem with constant coefficients on the half plane £ = {(x, y, z) £ M 3 : x > 0}. By rescaling and rotating the 
coordinates if necessary, one can bring the spacetime metric at p in the following form (see Q for details), 

ds 2 \ = -dt 2 + (dx + P x dt) 2 + dy 2 + dz 2 , 
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where /3 X represents the normal component of the shift vector at p. 2 As mentioned above, we set the normal component 
of the shift to zero since we arc interested in controlling the number of ingoing fields at the boundary. Therefore, we 
freeze the coefficients in front of the derivative operators in the evolution equations (JTHHJ) to the values 

a=l, f3 k =0, = 0, % = Sij. 

For the following, we use the standard operators from vector calculus grad , curl , div defined by 

(grad^ = d l( j), (curlX), = e m d k X l , div X = 8 k X k , 

for scalar and vector fields <f> and X, respectively. They satisfy the identities 

curlgrad0 = O, divcurlX = 0, div grad cf) = A(j>, (18) 
curl curl X — — AX + grad div X, (19) 

where A = d k d k denotes the standard Laplacian. We consider the following generalization from tensor calculus: 
(gradJQy := d (l X f) - ^S l0 d k X k , (curlT) tf := e kl{l d k T l jh (divT), := d%j, 

where T is a symmetric, traceless tensor field. Notice that by definition, grad X and curl T are symmetric, traceless 
tensor fields. The following identities generalize the previous ones from vector calculus: 

curl grad X = — grad curl X, (20) 

div curl T = -curl div T, (21) 

div grad X = ^AX + 1 grad div X, (22) 
3 

curl curl T = —AT + - grad div T. (23) 

Notice that the identities (|20I22[) imply that curl grad grad <fi = and div grad grad cj> = 2gradA$/3. With this 
notation, the evolution equations (JTHSJ) in the high-frequency limit are 

(24) 
(25) 
(26) 

1 4 \ 

+ -grad div p- -grad K J , (27) 

idivjS, (28) 
grad/3, (29) 

grad r 2grad grad <p — grad grad a, (30) 
4 

;rad div ft — — grad K, (31) 

o 

where /q, Go and Hq are the values of /, G and H frozen at the point p, and where for notational simplicity we omit 
the tildes over 7, A and T in what follows. The constraints in the high-frequency limit are 

H = idivT-4A0 = O, (32) 

M = div A- - grad if = 0, (33) 
3 

C = r-div 7 = 0. (34) 
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-foK, 


i\ = 


-Aa, 
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G B, 


B = 


H yA/3 


4> - 


-\ K+ 
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-2A + 2 


A = 


1 . 

-2 A7 + 


f = 


+ 



2 A redefinition of x is possible that makes this term vanish; however, this transformation does not leave the spacetime domain M = 
[0, T] X S invariant, but results in a domain with a "moving boundary". 
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The goal of this section is to construct boundary conditions at x = which ensure a well-posed Cauchy evolution for 
the above system, guaranteeing that the constraints propagate. Our construction is based on similar ideas described 
in [TH, HH and start with analyzing the propagation of the gauge-invariant quantities, namely the curvature variables 
consisting of the linearized Weyl tensor and the linearized constraint variables H, M. and C defined above. 



A. Propagation of the gauge-invariant quantities 



H 


= o, 


C 


= 2M, 


M 


= div£, 


£ 


= curlZ? 


B 


= +curl£. 



Here, we focus our attention on the evolution of the gauge-invariant quantities consisting of the constraint variables 
"H, Ai and C and the linearized Weyl curvature tensor. The latter can be divided in electric and magnetic parts, 
which in terms of the BSSN variables, read 

£ = A + grad grad a = — -A7 + grad T — 2grad grad 0, (35) 
B = curl A. (36) 

Using the identities (|20If2"3"f it is not difficult to see that the evolution equations (|25I28H3"T1 imply the following FOSH 
system for these gauge-invariant quantities, 

(37) 
(38) 
(39) 

^gradM, (40) 

(41) 

The following observation is important for the argument below: the gauge- invariant quantities H, Ai, C, £ and B 
are not independent of each other. Indeed, one has 2div£> = 2divcurlA = curldivA = curlAI, resulting in the 
"super-constraint" P := 2divS-curl7W = 0. Also, one has 2div£ = -Adiv7-8grad A0/3 + Ar + grad div T/3 = 
AC — 2grad"H/3, resulting in the super-constraint Q := 2div£ — AC — 2grad"H/3 = 0. Therefore, one can modify 
the system above by adding multiples of the constraint P — and Q = to the equations. 

As an example, we may use either one of the super-constraints P = or Q — in order to derive the following 
decoupled system for the constraint variables H, A4 and C: 

H = 0, (42) 
C = 2M, (43) 
M = AM. (44) 

Our boundary conditions must ensure that these constraints propagate, that is, initial data satisfying H = 0, M. = 0, 
C = must lead to solutions satisfying these constraints everywhere on the spacetime domain. A simple condition 
which implies this property is the imposition of a homogeneous Dirichlet condition on the momentum constraint 
variable Ai, 

M = 0, (45) 

where here and in the following, the symbol = refers to equality at the boundary surface x = 0. In the next section, 
we show that this condition also leads to constraint-propagation for the full nonlinear BSSN equations. 

Now we go back to the full propagation system for the gauge-invariant quantities (|37H41|) . and try to specify 
boundary conditions for this system which incorporates the three constraint-preserving boundary conditions (I45p . 
In order to determine which boundary conditions may be specified for this system, we analyze the corresponding 
principal symbol and its characteristics (see the appendix). The non-trivial part of the principal symbol is 





n(g>M , (46) 
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where n is a one-form and where we use the notation (n A B)ij := Ski(in k B l j\, (n <E> M)ij :— nuAdfi — ^Sijn k Mk, 
(n ■ £)j := n l Eij. Decomposing 

M = M\\ n + M±, 
3 

£ = -£|m n ® n + 2n ® + 
3 

B = -B\\\\n<»n + 2n®B\\± + B±±, 

into pieces parallel and orthogonal to n, the eigenvalue problem A(n)u = Xu yields the three scalar equations 

XM n =S m , X£ m =M lb AS|m=0, (47) 

the three vector equations 

13 1 

AM_l = £||_l, X£\\ ± = --nAB n _ L + -Mx, AB||_l = - n A £ {{x , (48) 

and the two tensor equations 

\£±± = -n A B ±± , XB ±± = n A £ ±± . (49) 
From this we obtain the following characteristic fields with zero speeds 

B|HI, Zf = B\^-\nhMx, 



and the ones with nonzero speeds, 



v (±1) = %±*l||. 

V ± ±1} = 4£ ll± T2nAB ll± ±3M ± , 

V ± f } = £±xTnAB X x, 

with corresponding characteristic speeds A indicated by the superscripts (±1) and (0). Maximal dissipative boundary 
conditions have the form of a linear coupling between the in- and outgoing fields corresponding to the outward unit 
one-form normal to the boundary, i.e. n = —d x (see the appendix). In our case, this means that we may specify 
conditions of the following form, 

= C V(-V + G, V[ +1) = Cl V ( - 1] + G X , V[ + V = C2 v[-V + GxXj (50) 

with coupling constants Co, c\ and C2 and boundary data G, Gj_ and G±±. In order to fix the coupling constants Co, 
ci and C2 we first remark that the fields Vj^ 1 ^ are related to the linearized Weyl scalars 

#o = C a p jS l a m^Pm s , * 4 = C a p 1 sk a fh k" ) fh 5 

in the following way: the Newman-Penrose null tetrad {l a , k a , m Q , to q } is defined by the time-evolution vector field 
d t and the normal to the boundary n = —d x according to 

1= — (dt + n), k = -j={d t -n), 

up to a rotation m — > e lv m. Then, we have, 

V^ 1 ^ = (ty rh ®rh + ^ m ® to), Wx = (^4 m ® m + ^4fh ® fri). 
Therefore, choosing C2 = in Eq. (|50[) we obtain the boundary condition 

Vl + V = G ±± , (51) 

which is equivalent to specifying the Weyl scalar ^0 at the boundary. In particular, one might choose 

Gj_j_ = V. I" which freezes <J/o to its initial value. This condition has been shown to yield a reflection coef- 

t=o 

ficient that decays as fast as (fci?)~ 4 for monochromatic gravitational radiation with wave number k and a spherical 
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outer boundary of radius R 0, @] ■ It has also been tested numerically @ for gravitational waves propagating about 
a Schwarzschild black hole and shown to outperform other currently used boundary conditions. 

Next, choosing cq = 1 and G = in Eq. (f50| . is equivalent to setting the constraint .Mm to zero at the boundary. 

On the other hand, the form of V[£ does not allow one to set the orthogonal component, M±, of the momentum 
constraint variable to zero at the boundary. 3 For this reason, we perform a slight modification to the propagation 
system (jBTHITj) by using the super-constraint P — 2divS — curlA-1 = 0. Namely, we replace Eq. (|^0|) by 



£ = — curlB + — grad M + n <g) 



a | dWB — -curl.M 



(52) 



where n = — d x is the outward unit one-form normal to the boundary. With this modification, the only change in 
the eigenvalue problem A(n)u = u is the middle vector equation in Eq. (|48l) which now simplifies to A£m_l = M±. 

Therefore, the characteristic fields remain unchanged except for the fields V r _[ ±1 ' ) which have to be replaced by 



V[ ±1) =E U ±M± 



(53) 



With this modification, it is now possible to impose the momentum constraint at the boundary by choosing cq = c\ = 1 
and G = and G± = in Eq. (I5TH) . However, it remains to prove that these boundary conditions, together with the 
condition (fSTj) on the Weyl scalar \&o are maximal dissipative and that the modified evolution system (|37I38I39I52[ 
@J} is still symmetric hyperbolic. For this, it is convenient to replace B by the new variable JC defined by 

JC :=B-n®(nhM). 

In order to write down the principal symbol of the resulting evolution equations, we choose standard Cartesian 
coordinates x,y,z on S such that n = —d x . The principal symbol with respect to an arbitrary one-form m = 
m x dx + mA<ix A then reads 



A(m) 



f &xx \ 

Mx 

£xB 
K-xB 

Mb 
£ab 

V &AB / 



I - £ AB mA fc Bx + mx M x - m A M A 
+£ AB m A £ Bx 

+ m £a X 

-£ CD m c K,DB - \e B C m c K, xx + m x M.B + \m B M x 
+e CD m c £DB + \e B c m c £xx 
m x £ xB + m A £ AB - \m B £ xx 



-m x £ C{A K. '° B) +e C ( A m c K. 
\ +m x e ciA £ c B) - £ C ( A m c £ B) 



B)x ~ \8 AB £ CD m c ]C Dx + m( A M B ) - \8 AB m c M c 

Cp ~ ^U AB £ CD m C £Dx 



(54) 



where the indices A, B, C, D refer to the coordinates y and z, and where we have defined £ AB := £ AB — hS AB 8 £q d 
and analogous for IC AB . It is straightforward to verify that this symbol is symmetric with respect to the symmetrizer 
H defined by 



U T HU = BL.+K2 



Mi 



28 AB £ xA £ xB + 28 AU K xA K xB + 2S AB M A M B 



+ 28 AlJ 8 BU £ AB £ CD + 28 au 8^!C A bICcd, 

where U — {£ xx , IC XX , M x ,£ x b, l^xB, Mb,£ A b,I^ A b) T , that is, HA(m) = A(m) T H for all one-forms m. In a 
coordinate-independent notation, this reads 



(7 T H[/=|£||||| 2 
In particular, we have, 

U T HA(n)U 



|/C, 



mi 



2|^||±| 2 



2|/C|| ± | 2 



-2£ XX M X - 4£ xA M A 
I (|y(+ 1 )| 2 - \v { - 1] \ 2 



A£ AB e AC t c l 



V 



(+1)|2 



-w 



-2\M_ 



(-D |2 



2\£ 



±±\ 



2\K 



(I 



3 The boundary condition M± = is equivalent to = 1 ' — 4n A Z± \ and involves a zero speed field. 

4 This modification is partially motivated by the "boundary adapted system" in Ref. Q. 
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For our choice of boundary conditions, = y( x ) , V^ +1 ) = V_[ ^ an( l V^ 1 ^ = G_l_l we obtain 

U T HA(n)U = |G ± J 2 - |y[- 1} | 2 < |G ± J 2 . 

From this, we see that our boundary conditions are maximal dissipative. 
We may summarize the results of this subsection in the following 

Lemma 1 The IBVP on the spacetime domain M := [0,oo)xE consisting of the evolution equations I{37\38139152\41\ ), 
initial data for T-L, C, Ai, £ and B in L 2 (E) and the boundary conditions 

M = Q, £±± — ra A B±± = G±±, 

where n = —d x is the unit outward vector to 5E and G±± lies in L 2 (T), T := [0, oo) x <9E, is well posed. 

The first boundary condition guarantees constraint-preservation, and the second condition specifies data to the 
Weyl scalar fy . In particular, we obtain an L 2 estimate for the gauge-invariant variables H, C, A4, £ and B. In 
view of Eqs. (|33l35l36p this yields L 2 estimates for A, curl A and div A provided we have appropriate estimates for 
the lapse, a, and the trace of the extrinsic curvature, K (see the next subsection). We might ask whether or not 
the L 2 bounds for curl A and div A arc sufficient to bound the L 2 norm of the full gradient of A. On K 3 , it can be 
proven that an L 2 bound on curl A and div A implies an L 2 bound on grad A However, this is not true in general 
on the half-plane E as the following example shows: Let \ be an arbitrary harmonic function on E which decays 
exponentially to zero as |x| — > oo and let A :— gradgradx- Then, curlyl = and div A = 2gradAx/3 = 0, but 
A^O unless boundary conditions at x = force \ t° be linear. 



B. Propagation of lapse and shift 



Next, we analyze the gauge sector, that is, the propagation of lapse and shift described in the small amplitude, 
high-frequency limit by the evolution Eqs. ([24H27]) . It is useful to introduce the quantities a := — / _1 grado:, 
D := Gq Miv (3, R := G " 1 curl/3 in terms of which these equations yield the following first-order system 

(55) 
(56) 
(57) 
(58) 

R + grad D — — grad 

Go 



a 


= -f K, 


$ 


= GqB, 


K 


= iodiva, 


a 


= grad K, 


is 


= (-~c 


R 


= curl_B, 


D 


= divS, 



(59) 

(60) 
(61) 

where kq := AGqHq/3 and where we have used the identity (|19p . As in the previous subsection, this system is subject 
to super-constraints, which are 

C a ■= a + / _1 grad a = 0, 
C D := D-Go 1 div^ = 0, 
C R := i?-G " 1 curl/3 = 0. 

As we show now, this system can be cast into FOSH form provided the condition 

«o ^ fo (62) 
holds. In order to see this, we rewrite the last three equations in terms of the new fields 

Go(jo-Ko) Go(/o-ko) 



and use the constraint C a = 0, obtaining 



C = K f-|Curli? + gradFj , (63) 

R = curlG, (64) 
F = divG. (65) 
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The system (|55I56I57I58I63I64I65P is symmetric hyperbolic; a symmetrizcr H for the principal symbol A(n) is given 
by the quadratic form 

v T Hv = a 2 + \(3\ 2 + K 2 + f a \a\ 2 + \C\ 2 + ^k \R\ 2 + k F 2 , 
where v — (a, /3, K, a, C, R, F) T . In terms of the characteristic fields 

.an 



y<±^:=A-±V9ba||, wf^:=C n ±^F, W± ±VU ^ ) := C± T \^ n A R, 



we have 



v T HA(n)v = 2f Ka\\ +2k q FC\\ + -k (C A R)\\ 



To 



|y(+v75)|2 _ Iy(-V75)|5 



\ W {+V*Z)\2_\ W (-V*0\2 



Therefore, we are allowed to impose the following boundary conditions at x = 0, 



ii 



(66) 



with coupling constants ao, ai and 02 smaller than or equal to one in magnitude and boundary data go, go and g±, 
where we recall that || and _L refer to components parallel and orthogonal to n = —d x , respectively. In terms of lapse 
and shift, this is equivalent to 



(oq - 1)q = \ffo(a + 1)91 1 a + fog, 



(01 - 1) P\\ ~ 



(« 2 - 1) A 



Ko 



fo - ko 



-dn 



K {ai + 1) div/3 



no a 
fa — «o fo 



G 



o9\\, 



/o — K a 



d±a = 



/3ko~ 



(a2 + l)(d\\Px-d ± p\\)-G g ± . 



(67) 
(68) 

(69) 



As mentioned in the introduction, we fix the normal component, /3n, of the shift to zero. This requires that ao = 
— a\ = 1 and g = g\\ = 0. Choosing also 02 = 0, 5 we obtain the following boundary conditions in the gauge sector, 



d\\a = 0, 



011=0, 



«0 



fo - «0 



d±a + G g±- 



(70) 



Summarizing, we have 



Lemma 2 The IB VP consisting of the evolution equations i f 551 56\57\ 58\ 6S\ 64] 65\) , initial data for a, f3, K , a, C , R 

and F in L 2 (T,) and the boundary conditions |7fl| ) with boundary data g± £Lr(T) is well posed. 

Furthermore, the solution satisfies the constraints grada = —foa, div/3 = GqF — KoK/(fo — K o), curl/3 = GqR if 
the initial data satisfy these constraints. 

Proof. The first part follows from the considerations above, since we have a FOSH system with maximal dissipative 
boundary conditions. For the second part, we notice that the FOSH system implies the following evolution equations 
for the constraint variables C a , Co and Cr: 



C a = 0, C 



D 



0, 



Cr 



nofo 



Go(fo - /to) 



curlC a . 



□ 



Here, different choices for are also possible. The choice ai = is motivated by requiring a Sommerfeld-like boundary condition for 
the orthogonal part, of the shift vector. A natural alternative is 0,2 = 1 which (together with the condition fin =0) controls the 
tangential components of the shift at the boundary. 
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In particular, we obtain I? estimates for lapse and shift and their derivatives a = —foK, $ = GqB, grad a — —foa, 
curl/3 = GqR, div /3 = GqD. With the boundary condition flu =0 this also implies an L? estimate for the symmetric 
and traceless gradient, grad (3, of the shift. This follows from the integral identity 

J |grad/3| 2 d 3 a; = j Q|curl/3| 2 + ||div/3| 2 ^ d 3 x - 2 J [3 n drv f3 x dydz 

£ £ 9£ 

which can be proven using integration by parts. Therefore, we obtain L 2 estimates for lapse and shift and all their 
first-order derivatives. Their second derivatives can be estimated by first taking time- and tangential derivatives of 
the evolution equations and boundary conditions and repeating the above analysis to obtain an L 2 bound for a and 
/3, d±a, d±/3 and their first derivatives and then estimating the second normal derivatives d 2 a and <9 2 /3 using the 

equations a — f Aa and (3 = kq(3A/3 + grad div (3 — 4grad K) /4. 



C. Reconstructing the metric and connection variables 

Finally, we discuss how the metric variables 7, <f> and the connection variables A and T may be reconstructed from 
the gauge- invariant quantities H, C, M., £ and B and lapse and shift. First, </> and T may be obtained after integrating 
the linearized BSSN Eqs. (|28|l and (j3"Tj) . respectively. Then, the trace- free part of the extrinsic curvature is obtained 
by integrating the equation A = £ — grad grad a in time. Then, 7 is obtained after integrating Eq. (f2T)f in time. It 
remains to show that the resulting fields yield a solution to the linearized BSSN system equations with the required 
boundary conditions. We state our main result in 

Theorem 1 (Well-posedness for the linearized IBVP) Let £ = {(x,y,z) G R 3 : x > 0} be the half plane, 
M := [0, 00) x S spacetime and T '■= [0, 00) x <9E the time-like boundary. Denote by n = —d x the unit outward vector 
to 9E. Let kq — 4GqHq/3 and suppose the condition k,q 7^ /o holds. 

Given initial data for the fields a, K , j3, B, <f>, 7, A, T in C^°(E) at t — 0, and given boundary data G±± and g± 
in Cq°(T), there exists a unique smooth solution of the linearized BSSN system ^24[\31\) on M satisfying the boundary 
conditions 

2 

div A grad if = 0, (71) 




Gxx, (72) 
0, (73) 
0, (74) 



/o - 



dj_a + G .g_L, (75) 



where £ = —^Aj + gradT — 2gradgrad</> and B — curl A are the electric and magnetic part, respectively, of the 
linearized Weyl tensor. 

Furthermore, if the initial data satisfy the linearized constraints % = 0, M = and C = defined in L32i33\3^\ ), 
then the solution automatically satisfies these constraints on M . 

Proof. Existence: In a first step, we define the gauge-invariant quantities %, Ai, C, £, B at t = by using the 
Eqs. (|32I33I34I35I36I) . respectively. By the assumption, they lie in C^°(E). These quantities can be extended to M 
by solving the well-posed IBVP described in Lemma [1] 

In a next step, we define the quantities a := — / _1 grad a, C := B + Gq 1 koo/(/o — Ko), F :— Gq 1 [div (3 + kqK/ (fa — 
kq)] and R := G ( ^" 1 curl/3 at t = 0, and extend the fields a, f3, K, a, C, R, F to M by solving the well-posed IBVP 
described in Lemma [2] Define B(t) := C(t) — Gq 1 Koa(i)/(/o — kq) for all t > 0. According to the statement of this 
lemma, the above definitions for a, C, F and R are preserved, i.e. they are valid on M. 
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Next, define for t > 0, 



r(t) 
lit) 



r(o) 



- g(Q) 

#0 



A(0) + y [£ (s) — grad grad a(s)] ds, 
o 

t 

0(0) -i / [*(*)- div0(*)]tfe, 



t 

7(0) -2 y L4(s)- grades)] ds. 



(76) 
(77) 

(78) 

(79) 



We claim that the smooth fields a(t), K(t), f3(t), B(t), <j>(t), 7(i), A(t) and T(t) obtained in this way solve the linearized 
BSSN equations ([24H3T]) . For this, we first claim that idivT 4A</> = U, div A - fgradif = M, T — div7 = C, 
— 5A7 + gradr — 2gradgrad</> = £ and that curl A — B for all t > 0. In order to prove this, we first notice that 
these equalities are true for t — by definition of the gauge-invariant quantities %, Ai, C, £, B at t = 0. Then, the 
Eqs. (|37I38I39I41I52I57I59|) and the above definitions for T, A, <f> and 7 imply that 



dt 

^ [M-divA + |grad K ) = 0, 

<xr \ 3 / 

— (C - T + div7) = 2 |X - divA + -gradK 

<k \ 3 

^ £ + - A7 — grad r + 2grad grad < 

3 / 2 

= —curl (B — curl A) + -grad f M. — div A + -grad K 



■ n <g> 



n A divS curl M 

2 



dt 



(B - curl A) = 0, 



d ( 1 
— divS curl.M 

dt\ 2 



0. 



(80) 
(81) 
(82) 

(83) 
(84) 
(85) 



which proves the first claim. Now it is straightforward to verify that the fields a(t), K(t), (3(t), B(t), cf>(t), 7(t), A(t) 
and r(t) satisfy the linearized BSSN equations (l24T[3"Tj) . 

Uniqueness: Suppose flW, </>W, 7 «, r« and a< 2 ), , /?( 2 \ B®, </>( 2 ), 7 ( 2 ), A< 2 ), T< 2 ) 

are two smooth solutions of the linearized BSSN equations (|2"4"H3"Tj) which are identical on the initial surface t = 
and which both satisfy the boundary conditions (ITTHTS")) with identical data Gj_j_ and g±. Then, a := — 
(3 := /3' 2 ' — fi( x \ ... , r := r' 2 ) — r' 1 ) is also a solution with trivial initial data satisfying the boundary conditions 
(|71H75|) with homogeneous boundary data. In particular, the associated gauge-invariant quantities %, Ai, C, £, B 
defined by Eqs. (|32I33I34I35I36P satisfy the IBVP described in Lemma[T]with trivial initial and boundary data. Since 
this IBVP is well posed, it follows that these quantities are zero on M. Similarly, it follows from the well-posedness 
of the IBVP described in Lemma [2] that a, f3, K and B must be zero on M. Next, it follows from Eqs. ([281311) that 
and r are zero on M. Next, since = -2£ = A7, it follows from Eq. that A = on M. Finally, Eq. (|2"9"|) 
implies that also 7 = on M which shows that the two solutions are identical. 

Constraint preservation: This will be proven in the next section for the more general case of the full nonlinear 
BSSN equations. rj 



D. Geometric uniqueness 



In this subsection we analyze the residual gauge freedom in our initial-boundary value formulation of the linearized 
BSSN system. With respect to an infinitesimal coordinate transformation generated by a vector field X^, say, the 
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linearized four metric h^ v transforms according to 

hftv i-> h'^ = + £ X V^ = + d^X v + d v X^, 

where is the Minkowski metric. This induces the following transformations on the linearized BSSN variables, 

a' = a + F, K' = K - AE, /?' = /? + Y, B' = B + Z, (86) 

<t>' = <j>+-divX, 7' = 7 + 2gradX, A! = A — grad grad E, (87) 
6 

r' = T + AX + igrad div X, (88) 
3 

where we have defined E := X°, X := (X\X 2 ,X 3 ) and F := E, Y := X - grad£, Z := G^Y . The linearized 
BSSN equations (f24l f3T]) are invariant with respect to these transformations if and only if E, F, Y and Z satisfy the 
equations 

E = F, (89) 

F = f Q AE, (90) 

Y = G a Z, (91) 

Z = H [AY + igrad div Y + -^grad F) . (92) 
V 3 3/o/ 

Similarly, the boundary conditions (|71H75|) are invariant provided that F and Y satisfy 

d l{ F = 0, Yj|=0, (d t + ^-a\^Y ± = - r ^—a ± F + G w ± , (93) 

V ^ / Jo - k 

where w± = g' ± — g± and g± and g' ± are the boundary data corresponding to the two different gauges. Since the 
Eqs. dBHHHl]) yield the two second-order equations 

F = foAF, Y = G H (AY + igrad div Y + ^Srad F^j (94) 

which are identical to the second-order equations for a and /3 obtained from Eqs. (f2"4Tf2"?j) . and since the boundary 
conditions are identical to the boundary conditions (I73"HT5)) , it follows from the results in Lemma[5]that the IB VP 
for F and Y consisting of the evolution equations (l94l) . the boundary conditions (|93|) . initial data for F, F = foAE, Y, 
Y = GqZ in L 2 (Yj) and boundary data w± in L 2 (T), is well posed. Therefore, the residual gauge freedom is entirely 
determined by the initial values for E, F, Y, Z and the boundary data g± and g' ± . This allows us to formulate the 
following result: 

Theorem 2 (Geometric Uniqueness) Two smooth solutions of the IB VP described in Theorem[J\ are related to 
each other by an infinitesimal coordinate transformation if and only if there exists two smooth functions F and E on 
£, three smooth vector fields X , Y and Z onE and a smooth vector field w± on T such that the initial data of the two 
solutions are related to each other by the transformations i86\87]8W and the boundary data (G_u_, gx)> (Gj_x>9±.) * s 
related to each other by 

G'±j_ = Gj_j_, g' ± = gx + wx- 

Proof. The "only if part of the proof follows from the considerations above. For the "if part, we use the fields E, 
F, Y and Z to set up initial data for F, F = foAE, Y, Y — GoZ, and the field wx to set up boundary data for the 
IB VP described by Eqs. (IM|) and (li?3"|) . Since this problem is well posed, we may extend the fields F and Y uniquely 
to all M. Defining 

t t 
E{t) := E(0) + J F(s)ds, X(t) := X(0) + j [Y(s) + grades)] ds, 


we obtain a vector field (X^^t)) = (E(t),X(t)) on M which, by construction, generates an infinitesimal coordinate 
transformation between two solutions of the linearized IBVP for the BSSN system. rj 
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IV. CONSTRAINT PRESERVATION 



In this section we prove that smooth enough solutions satisfying the nonlinear BSSN equations <jT] [8]) and our 
boundary conditions satisfy the nonlinear constraints (|11I12I13[) everywhere on the spacetime domain M — [0,T] x Y, 
if the initial data satisfy these constraints on the initial slice {0} x E. The proof is based on the FOSH system derived 
in [29|, H3| governing the propagation of the constraint fields with the parameter choice m = l, a = 1/2, The detailed 
equations describing this system will not be needed for the following argument. Only its form, 



d C 



a 



[A(uYd l C + B(u)C] 



(95) 



is important. Here, C — (C l , %, S := H + Z/2, Mj, Zuj\, Z^) are the constraint variables, where we have decomposed 
the variable := (diC k )~fkj — Z(ij) + +JijZ/3 into its trace-free symmetric part, Zuj\, its antisymmetric part, 
Zuj], and its trace, Z = r ) tJ Zij, where u = [a, ft 1 , <p, K,^ij, Aij) are the main variables, and where A 4 , i = 1,2,3, and 
B are matrix-valued functions of u. The principal symbol A(n) = K(u) l rii is given by 



A(n) 



H 

S 

\ %]) / 



/ \ 



v? Mj 

\ nj S + \n l Z {l2) + WZ M 
2(n (i M i) ) Ti ^ 



(96) 



Here n l = ~f lJ rij and ni is normalized such that n^n 1 = 1. This system is symmetric hyperbolic; a symmetrizer 
H = H T is given by the quadratic form 

C T HC = +n 2 + ^S 2 + fiMiMj + 1 7 <fc 7 ^Z (ij) Z (H) + ~ 7 ik 7 jl Z m Z [kl] . 

In order to obtain the class of maximal dissipative boundary conditions, we find, following the appendix, the charac- 
teristic speeds and fields which are defined as the eigenvalues and corresponding projections of C onto the eigenspaces 
of A(n). They are given by 



n = o, c k 



/' 



n, 



s 



i?Z, 



±1, V {±1) = n j Mj ± 



\n l n 3 Z t 



( Z (kj) 



M = ±l, W. 



ik] 



n 



Z[ki]^) 



Z(ki)P kl ij , 



where III and u are the projection operators defined in Table HI In terms of these fields we find 



C J HA(n)C = - {V i+1) f - {V^f 



which shows that the boundary conditions 



-7 

2 ' 



1 J * J 



w. 



(+1) 



b 2 W, 



(-1) 



(97) 



where b\ and 62 are two functions on the boundary T whose magnitude is smaller than or equal to one, and where 
rii refers to the outward unit normal to 9S, are in maximal dissipative form. The particular case b\ = 62 = — 1 
corresponds to imposing the momentum constraint on the boundary, M.j = 0. Since the constraint evolution system 
(|95p is linear and homogeneous, this implies that the unique solution with trivial initial data is zero. Therefore, 
imposing the momentum constraint on the boundary guarantees that solutions which satisfy the constraints initially 
automatically satisfy them on M. Furthermore, the energy estimate from the appendix shows that a suitable norm 
of the constraint violations which are usually present in numerical calculations is bounded by a (time-dependent) 
constant times the norm of the initial constraint violations. Therefore, for each fixed time t, smaller and smaller initial 
constraint violations lead to smaller and smaller constraint violations at t. Finally, we mention that our boundary 
conditions do not modify the key property of the turducken approach (29l . [30l | regarding the causal propagation of 
constraint violations. In particular, constraint violations initiating inside a black hole cannot propagate outside the 
hole for the parameter choice m = 1 used in the present work. 
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V. A NOTE ON THE NUMERICAL IMPLEMENTATION 

Regarding the numerical implementation of our boundary conditions we have the following comments: from the 
system of evolution equations, Eqs. ([THE]), it is clear that a question arises at the moment of discretizing the spatial 
derivatives at the boundary points. For finite difference methods, one possibility is to use one-sided differences at 
those points. However, this simple method cannot yield a convergent numerical schemes in general since it does not 
incorporate the information about the boundary conditions, and so the corresponding system at the continuum is 
underdetermined. (Using one-sided differences would be enough if the boundary is purely outflow, as happens at 
the inner boundary of some excision schemes.) The question then is how to incorporate the boundary conditions 
summarized in Table [X] into the finite difference discretization. Without specifying a detailed method for doing so, let 
us nevertheless mention a few ideas. 

Concerning the gauge boundary conditions d n a = 0, j3 n — on the lapse and the normal component of the shift, 
they could be used in order to extrapolate the values of a and (3 n to ghost zones outside the numerical grid. These 
values could then be used in order to compute centered differences at the boundary points. A similar procedure could 
be applied to the gauge condition for the tangential components of the shift. 

In order to implement the constraint preserving and specifying conditions, first note that they are somewhat 
less common than the gauge boundary ones, as they involve second-order derivatives of the metric fields. In this case, 
we propose to follow a procedure described in |35|, where constraint preserving boundary conditions, including the 
specifying condition, are numerically implemented and tested for a first-order symmetric hyperbolic formulation 
of Einstein's equations. Applied to the BSSN system, this method consists in adding terms to the right hand side 
of the evolution equations, Eqs. ([THE]): which are proportional to the expressions 

SirGNji and = P v i m £ij — (n l P 3 i m — n P 13 im) &kij — P l3 imGij> defining the momentum constraint and the 
$0 specifying boundary conditions, respectively The proportionality factors are then chosen in such a way that 
normal derivatives corresponding to the combinations of the characteristic fields appearing in the expressions for $^ 
and ^ij, are eliminated. The characteristic fields for the BSSN Eqs. ([THE]) have been calculated in [1J] based on the 
pseudo-differential approach described in [3ll . |32| . The ones that are relevant to our problem are 

where n refers to the unit outward normal to the boundary, and the projectors IT*- and P 13 ' m are defined in Table [U 
In terms of these fields, we have 

*j = \d n (> 1} + v^) n 3 + \d n + v^) + (98) 

= M^ X) + ■■•> (") 

where ... refers to terms containing tangential derivatives or undifferentiated combinations of the characteristic fields 
only. Therefore, at the boundary points, $j and have to be added to the evolution equations in such a way that 

the normal derivatives of + w^ 1 - 1 , Vj +1 ^ + v 1 ^ 1 ^ and v^j are eliminated. 

Once this calculation has been performed and the coefficients in front of and have been determined, the BSSN 
evolution equations ([THE]) with the additional boundary terms may be discretized using finite difference operators. At 
the boundary, one-sided differences might be used for all fields except for lapse and shift for which centered differences 
with the extrapolated values at the ghost zones are implemented, as described above. We stress that once the above 
mentioned coefficients have been determined, the characteristic fields are not required anymore and the system might 
be implemented numerically. 

For a different method for discretizing constraint-preserving boundary conditions based on pseudo-spectral methods 
we refer the reader to 1361. 



VI. CONCLUSIONS 

In this work we have analyzed the IBVP for the BSSN evolution system of Einstein's field equations with a 
hyperbolic X-driver and Gamma-driver condition [17[ for lapse and shift, as used in current numerical simulations. 
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Unlike the harmonic formulation which has been motivated by the mathematical structure of the equations and the 
understanding of the Cauchy formulation in General Relativity, the BSSN system has been developed and improved 
based on its capability of numerically evolving binary black hole spacetimes in a stable way. Therefore, it is not 
evident at all that mathematical questions like the well-posedness of the IBVP may be answered. 

One of the important results obtained in this article is that the BSSN evolution system yields a "nice" evolution 
system for the constraint variables and the Weyl curvature and a "nice" evolution system for lapse and shift. These 
two systems are "nice" in the sense that they may be cast into first-order symmetric hyperbolic form, from which 
one concludes that they yield a well-posed time evolution with standard prescriptions for constructing boundary 
conditions. At the linearized level, the two system decouple from each other and the gauge-dependent fields may be 
constructed from their solutions as we have shown in Theorem [1] At the nonlinear level, these systems couple to 
each other through the metric fields, and one needs to consider the coupled system together with the evolution of the 
metric fields. If this large system could be cast into first-order symmetric hyperbolic form as well, a well-posedness 
proof would follow from the standard theorems quoted in the appendix. We have not investigated this issue in the 
present article. Instead, we have extrapolated the boundary conditions constructed in the linearized case to the full 
nonlinear BSSN system. These conditions are summarized in Table HI and they are conjectured to yield a well-posed 
system in the nonlinear case as well. 

Summarizing, our new boundary conditions for the nonlinear BSSN system have the following properties: they 
preserve the constraints throughout evolution, as shown in Sec. IIV1 they yield a well-posed IBVP in the linearized 
case, and furthermore, they control the Weyl scalar a t the boundary, a condition that is currently used in binary 
black hole simulations based on the harmonic formulation of Einstein's equations (see, for instance, Ref. [371 ]) and has 
been tested for its accuracy for nonlinear waves on a Schwarzschild background [9j . 

In a first step towards their actual numerical implementation, we have sketched a possible numerical method for 
imposing our boundary conditions which is based on finite-differences and projection techniques. It is the hope of 
the authors that such an implementation eventually improves the accuracy of the simulations and avoids the need 
of pushing the boundary so far away that it is causally disconnected from the region where physics is extracted. 
Our boundary conditions may also be useful for matching the solution at interface boundaries, or for the Cauchy- 
characteristic matching approach. 
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In this appendix we summarize some known results about FOSH systems with maximal dissipative boundary 
conditions [38|, [3!|. Let £ C R™ be an open subset of M™ with closure E and C°° boundary <9£. We consider the 
following IBVP on the spacetime region M := [0, oo) x E, 



Here, A 1 , A 2 , ... ,A n , D : M -> Mat(m x m, E) are C°° matrix-valued functions on M, F : M -> R m is a C°° vector- 
valued function on M and B : T — > Mat(m x m, R) is a C°° matrix-valued function on the boundary T '■= [0, oo) x 9E. 
The data consists of the initial data uq E L 2 (T,) and the boundary data G E L 2 (T), and u : M —> R m is the solution 
vector, lying in an appropriate function space. 

We require the following conditions on the principal symbol A(t,x,k) :— A l (t,x)ki, k — (fci, &2, k n ) E R™, and 
the boundary matrix B(t,x): 

(a) There exists a C°° matrix-valued function H : M — > Mat(m x m,R), called the symmetrizer, such that 

(i) H(t,x) = H(t,x) T is symmetric for all (t,x) E M. 

(ii) There exists a constant C > such that C~ l \u\ 2 < u T H(t, x)u < C\u\ 2 for all (i, x) E M and all u E W n , 
that is, H is uniformly positive definite. 
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Appendix A: Symmetric hyperbolic first-order systems with boundaries 



d t u = A i (t,x)d i u + D(t,x)u + F(t,x), (>0,ieE, 
u(0,x) = Uo(x), x E E, 

B(t,x)u(t,x) = G(t,x), t>0,xe<9E. 



(Al) 
(A2) 
(A3) 
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(iii) H(t, x)A(t, x, k) = A(t, x, k) T H(t, x) is symmetric for all (t, x) € M and all k € R™. 

(b) Let J(t, x) := H(t, x)A(t, x, k), where (t, i)eT and where k denotes the outward unit normal to 9E. Then, we 
assume that for each p = (t, x) € T the boundary space 

V p := {u £ K m : B(t, x)u = 0} C M m 

is maximal non-positive with respect to J(t,x), that is, 

(iv) u T J(t, x)u < for all u eV p . 

(v) V p is maximal with respect to the condition (iv), that is, if W D V p is a linear subspace of IR m containing 
V p which satisfies (iv), then W — V p . 

Furthermore, we require 

(vi) The rank of J(t,x) is constant. 

Under these conditions, it can be shown that the IB VP (|A1IA2IA3P is well posed [H-[4l[. In [H a generalization to 
other domains E with boundaries which might contain corners is also given. Generalizations to quasi-linear systems, 
where A l {t,x, u), D(t,x,u) and B(t,x,u) might also depend on the solution vector itself are also possible |42| . 

The meaning of the conditions (i-iv) is that they yield an a priori energy estimate. In order to see this, let u be a 
solution of the IB VP (|Al|A2|A3[) . By redefining u if necessary, we might assume that G = 0. Next, define the energy 
norm 



E{t) := J u(t,x) T H{t,x)u(t,x)d n x = \\H 1/2 (t,-)u(t,-)\\ 2 L2(s) , I > 0. 

E 



Taking a time derivative and using the symmetry of H and the evolution equations (|A1[) we obtain 

^E(t) = J ^2u(t,x) T H(t,x) [A i (t,x)d i u + D(t,x)u + F(t,x)] + u(t, xf[d t H{t, x)]u(t, x)}<Px. 
s 

Next, using the symmetry of H(t, x)A l (t, x) and Gauss's theorem, this yields 

^-E(t) = J u(t,x) T J(t,x)u(t,x)da + J u(t, x) T L(t, x)u(t, x)d n x + 2 J u(t,x) T H{t,x)F(t,x)d n x, 
as E E 

where da denotes the area element on 9E and where L(t,x) := dtH(t,x) — di[H(t,x)A l (t,x)} + H(t,x)D(t,x) + 
D(t,x) T H(t,x). The first term on the right-hand side is non-positive due to condition (iv). Using Schwarz' inequality 
and the condition (ii) the second and third term on the right-hand side can be estimated as 

u(t, x) T L(t, x)u(t, x)d n x < a(t)\\u(t, .)||| 2(s) < Ca(t)E(t), 

where a(t) :— sup{||L(s,x)|| :0<s<t,i€ £}, and 

2 J u(t,x) T H(t,x)F(t,x)d n x<2E(t)VC\\F(t,.)\\ L 2 {l:) < ^L E (t) + K 2 \\F(t, .)\\ 2 L2{ ^, 

E 

where K > is a constant, respectively. Summarizing, we obtain 

±E(t)<b(t)E(t)+K 2 \\F(t,.)\\ 2 L2{s) , 

where b(t) := C[a(t) + K~ 2 ]. Finally, using Gronwall's lemma we have 

t 

E{t) < e b(t)t E(0) + K 2 J e b{t){t - s) \\F{s, .)||| 2(s) ds, t > 0. (A4) 
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Such or similar a priori estimates are the basis for showing different properties of the solutions, like their uniqueness, 
continuous dependence on the data and the finite speed of propagation. The condition (v) is important for the 
existence of solutions, ensuring that not "too many" boundary conditions are specified. 

A more practical formulation of the conditions (iv) and (v) is the following: let p = (t,x) € T be fixed, and let k 
denote the outward unit normal to dT, at x. Since A(t, x, k) is symmetric with respect to the scalar product on R m 
defined by < u, v >:= u T H{t,x)v, u,v £ R m , there exists a basis ei,e2,...,e m of eigenvectors of A(t, x, k) which is 
orthonormal with respect to < •, • >. Let Ai, A 2 , X m be the corresponding eigenvalues, where we might assume that 
the first r of these eigenvalues are strictly negative, and the last s are strictly positive. We can expand any vector 
u € M. m as 

m 
3=1 

The coefficients Uj are called the characteristic fields; the associated speeds are given by Xj. Then, we have 

m r m 

u T J(t,x)u —< u, A(t,x,k)u >—^^\jU? = — ^^|Aj|u? + ^ 

j—1 j=l j—m — s+1 

In our work, the nonzero eigenvalues come in pairs, such that r = s and Ai = — A m , A2 = — X m -i, ■■■ >-V = — A m - r +i. 
In this case the above simplifies to 

u T J(t,x)u = \Xj\ (j^ (+) | 2 - l w j r) | 2 ) , 
3=1 

where we have defined := u m _j_i and Vj ' := uj. Clearly, the boundary conditions Vj + ^ = c^Vj ^ with 

|c^| < 1, j = 1,2,..., r, satisfy the conditions (iv) and (v). 
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